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ABSTRACT 

This paper considers the incorporation of constraints to enforce physically based conservation laws in 
the ensemble Kalman filter. In particular, constraints are used to ensure that the ensemble members 
and the ensemble mean conserve mass and remain nonnegative through measurement updates. In 
certain situations filtering algorithms such as the ensemble Kalman filter (EnKF) and ensemble 
transform Kalman filter (ETKF) yield updated ensembles that conserve mass but are negative, even 
though the actual states must be nonnegative. In such situations if negative values are set to zero, or 
a log transform is introduced, the total mass will not be conserved. In this study, mass and positivity are 
both preserved by formulating the filter update as a set of quadratic programming problems that in- 
corporate nonnegativity constraints. Simple numerical experiments indicate that this approach can 
have a significant positive impact on the posterior ensemble distribution, giving results that are more 
physically plausible both for individual ensemble members and for the ensemble mean. In two exam- 
ples, an update that includes a nonnegativity constraint is able to properly describe the transport of 
a sharp feature (e.g., a triangle or cone). A number of implementation questions still need to be 
addressed, particularly the need to develop a computationally efficient quadratic programming update 
for large ensemble. 


1. Introduction 

The importance of respecting physical conservation 
principles has long been recognized in numerical weather 
prediction modeling (Arakawa 1972; Arakawa and Lamb 
1977; Sadourny 1975; Janjic 1984; Janjic et al. 2011; Janjic 
and Gall 2012). One of the most basic of these principles 
is the need to conserve the total mass of air, water in 
its different phases, and relevant chemical species, in 
the latter two cases accounting properly for sources 
and sinks while maintaining the proper sign (positivity 
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or, more strictly, nonnegativity) in every grid volume 
of a computational model (e.g., Lin and Rood 1996, 
p. 2047). Smolarkiewicz and Margolin (1998, p. 460) 
write in a review paper that “the preservation of sign 
during numerical advection is the essential aspect of the 
stability and accuracy in modeling water phase -change or 
chemical processes.” Sign preservation should be an im- 
plicit requirement in any numerical algorithm that at- 
tempts to conserve mass. Although most of the concepts 
described in this paper apply to the conservation of other 
quantities, such as angular momentum and energy, we 
focus on mass conservation to make the discussion more 
specific and to illustrate ideas with simple examples. 

When ensemble Kalman filters are used for data as- 
similation two distinct mass conservation issues arise. 
Mass should be conserved in each member (replicate) of 
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the forecast ensemble produced by propagating states 
over time with a forecast model, and it should be con- 
served by the procedure used to update the ensemble 
members with observations. Conservation of mass dur- 
ing the forecast is dependent on the time- and space- 
discretized method used to obtain a numerical solution 
to the governing continuum equations. Although mass is 
conserved by construction in the original differential 
equations, it may not be conserved in the numerical 
solution. 

Similarly, conservation of mass during the ensemble 
filter update depends on the particular numerical method 
used to generate the replicates of the updated ensemble 
and to generate the mean of this ensemble (analysis 
itself). Mass may not be conserved during the update 
even when the numerical forecast technique is mass 
conservative. A number of methods have been proposed 
to deal with this issue. For example, Jacobs and Ngodock 
(2003) noted that mass in a simplified ID ocean model 
can be conserved when the model error in a representer 
algorithm is expressed in terms of the mass flux due to 
uncertainty in ocean depth rather than as additive error 
in the continuity equation. In land surface hydrology, 
Pan and Wood (2006) showed how to ensure conserva- 
tion of total water mass by imposing it as a “perfect 
observation” in a two-step Kalman filter approach. In an 
ocean data assimilation system, Brankart et al. (2003) 
imposed conservation of total mass, including positive 
layer thicknesses, through an a posteriori adjustment to 
the analyzed state. Positivity can also be ensured by 
introducing a change of state variables, using techniques 
such as Gaussian anamorphosis (e.g., Simon and Bertino 
2009). In atmospheric data assimilation, nonnegativity 
of the specific humidity has been imposed as a weak 
constraint in a three-dimensional variational data as- 
similation (3D-Var) implementation (Liu and Xue 2006; 
Liu et al. 2007). 

It is reasonable to ask if we should conserve mass in an 
ensemble Kalman filter update if the total mass in the 
system is uncertain. Measurements introduced during 
the update may provide useful information about this 
uncertain mass. If so, this information should be used to 
adjust and improve mass estimates. But this does not 
change the fact that a filtering algorithm should be able 
to preserve a known value of total mass through both the 
forecast and update steps. If the total mass is, in fact, 
unknown then it should be treated as uncertain and es- 
timated as part of an otherwise mass-conservative fil- 
tering procedure. Thus, we distinguish the need to 
respect conservation laws from the need to properly 
account for uncertainties in conserved quantities. In the 
simple examples considered here total mass is a known 
constant that does not need to be estimated, but this 


restriction could be relaxed by including uncertain 
sources/sinks for instance in the state vector so they can 
be updated when new information becomes available. 

Going beyond mass conservation, Cohn (2009) has 
shown in the context of minimum variance state esti- 
mation that conservation of total energy requires in- 
cluding a special term in the evolution equation for the 
state estimate that couples state and covariance evolu- 
tion. This requirement was subsequently formulated 
more generally as the principle of energetic consistency , 
and was used to study certain pathological behavior of 
ensemble-based data assimilation schemes (Cohn 2010). 
Cohn (2010) shows that the mild energy dissipation 
typical of numerical weather prediction models can lead 
to ensemble collapse through a feedback mechanism 
introduced by the assimilation of observations, and that 
this dissipative behavior can in principle be eliminated 
by including an appropriately scale-selective, anti- 
dissipative operator in the formulation of the ensemble 
data assimilation scheme. Using such an operator to 
maintain ensemble spread is a generalization of the co- 
variance inflation technique now commonly used in 
ensemble data assimilation schemes (Anderson and 
Anderson 1999, p. 2747). 

One objective of data assimilation is to use observa- 
tions to correct forecast errors in the vicinity of well- 
defined natural features such as fronts, filaments of 
chemical constituents, or plumes from surface emissions 
of aerosols. But data assimilation schemes have not 
traditionally been explicitly formulated to preserve such 
sharp features and, in fact, they often tend to blur and 
distort sharp interfaces (e.g., Lawson and Hansen 2005). 
Riishpjgaard (1998) proposed a type of state-dependent, 
anisotropic covariance modeling as a simple and direct 
analysis approach in the presence of sharp features; see 
Liu and Xue (2006) for an implementation. Hoffman 
et al. (1995) approached feature analysis by defining 
nontraditional, feature-based measures of spatial fore- 
cast error and minimizing them explicitly in variational 
data assimilation; see Gilleland et al. (2010) for a review 
of forecast verification methods utilizing such measures 
of forecast error. Lawson and Hansen (2005) showed 
that the performance of an ensemble Kalman filter in 
the presence of a well-defined feature can be improved 
dramatically by the use of alternative error models to 
redefine the state estimation problem. 

All of these studies highlight the importance of con- 
serving known mass, maintaining nonnegativity, and 
preserving feature geometry during data assimilation. In 
an ensemble context it can be argued that individual 
ensemble members as well as the ensemble mean should 
meet all of these requirements. Here we examine the 
issues of mass conservation, nonnegativity, and feature 
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preservation with two computational experiments that 
focus on simple features with well-defined shapes. In 
addition, we propose a new sequential ensemble data 
assimilation algorithm that enforces mass conservation 
and maintains nonnegativity by adding constraints to the 
ensemble Kalman filtering update. 

We can view the classical Kalman update at any given 
time as either an unconstrained regularized minimum 
variance estimator or, if we take a Bayesian perspective, 
as a method for generating the maximum a posteriori 
(MAP) estimate when the a posteriori probability is 
Gaussian. In either case, the update requires the solution 
of an unconstrained optimization problem. When the 
measurement operator is linear the problem objective 
function depends quadratically on the analysis (which is 
the decision variable to be determined). In this case the 
unconstrained optimization problem has a closed form 
solution given by the classical Kalman filter update 
equation. 

If the Kalman update does not satisfy physically based 
mass conservation or nonnegativity conditions, it is rea- 
sonable to enforce these conditions by adding appropri- 
ate constraints to the original unconstrained problem 
(Simon and Simon 2005). When the constraints are 
linear equalities and/or inequalities and the objective 
function Hessian is positive definite the resulting con- 
strained optimization problem is a convex quadratic 
program with a unique global minimum. It is important 
to emphasize that this constrained optimization problem 
is a quadratic program even if the forecast model is non- 
linear, so long as the constraints and observation operator 
are linear. 

It is possible to apply similar concepts to ensemble 
filtering problems, where we can obtain each member 
of the analysis ensemble by solving a replicate-specific 
quadratic programming problem, replacing the forecast 
mean with one of the forecast ensemble members and 
adding random measurement errors to the actual ob- 
servations. This procedure is analogous to the ensemble 
Kalman filter formulation proposed by Burgers et al. 
(1998). 

In our extension of the Burgers et al. (1998) ensemble 
Kalman update we include nonnegativity constraints 
when computing the replicates of the analysis ensemble. 
This ensures that the ensemble members are all physi- 
cally plausible, making it more likely that sample sta- 
tistics, such as the covariance of this ensemble, are also 
physically realistic. From a Bayesian perspective, the 
nonnegativity constraints provide additional prior in- 
formation that is not included in the classical formula- 
tion of the ensemble filtering problem. Consequently, 
the analysis replicates are no longer Gaussian but are 
strictly nonnegative. Our formulation of ensemble Kalman 


filtering as a sequence of static linear estimation problems 
makes the incorporation of physically based constraints 
a natural extension of the classical algorithm. 

It is useful to briefly distinguish the ensemble qua- 
dratic programming approach proposed here from other 
approaches that share some of its features. The review 
article of Simon (2010) gives an overview of various 
methods for incorporating constraints into the classical 
Kalman filter, including a version of quadratic pro- 
gramming. These classical methods are able to maintain 
mass conservation through the filter update but not 
during the forecast if the system dynamics are nonlinear. 
Our quadratic programming approach uses an ensemble 
forecast that is able to conserve mass and nonnegativity 
through the forecast step for each ensemble member, 
even for nonlinear problems, if the forecast model is 
properly formulated. 

Another ensemble data assimilation technique known 
as randomized maximum likelihood (RML) also solves 
an ensemble of optimization problems (Gu and Oliver 
2007; Emerick and Reynolds 2013). In this case each 
problem minimizes the batch mean-squared measure- 
ment misfit computed over all measurement times (per- 
haps with an additional quadratic regularization term) 
for a particular replicate. The forecast and observation 
models are formulated as nonlinear equality constraints 
and are incorporated into each optimization problem 
through a set of derived objective function gradients. 
The problem solution is obtained with an unconstrained 
nonlinear programming procedure. Our approach uses 
a different objective function at each measurement time 
as well as for each ensemble member since it is formu- 
lated as a sequence of static updates rather than as a batch 
algorithm. This makes it possible to formulate the opti- 
mization problem as an efficient quadratic programming 
problem that readily accommodates inequality con- 
straints. It is also compatible with the time-recursive 
structure of the classical ensemble Kalman filter. 

Other options for incorporating ensemble informa- 
tion include a number of different ensemble variational 
or hybrid data assimilation algorithms (see, e.g., Hamill 
and Snyder 2000; Lorenc 2003; Buehner 2005; Zupanski 
2005; Wang et al. 2007b, 2008; Wang 2010, 2011; Isaksen 
et al. 2010; Bonavita et al. 2012). These generally use 
a forecast ensemble to construct the hybrid covariance 
needed for a variational update. The decision variables 
in the optimization problem can include the analysis 
mean (Wang et al. 2008; Wang 2010, 2011) or individual 
ensemble members (Hamill and Snyder 2000; Isaksen 
et al. 2010; Bonavita et al. 2012). Additional inequality or 
equality constraints such as those used in our quadratic 
programming approach are typically not included in these 
hybrid methods. Despite these differences, our data 
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assimilation procedure shares important features with 
randomized maximum likelihood and ensemble varia- 
tional methods (Hamill and Snyder 2000; Zupanski 
2005) and it is possible to imagine variants that combine 
aspects of all three approaches. 

In this paper we compare our constrained Kalman 
filtering algorithm to an ensemble Kalman filter (EnKF; 
Evensen 2009; Houtekamer and Mitchell 1998; Burgers 
et al. 1998). Section 2 provides background and con- 
siders mass conservation and sign preservation for en- 
semble Kalman filters. Section 3 describes two numerical 
examples that illustrate the consequences of mass bal- 
ance errors in this class of filters, using results obtained 
with the EnKF. The first example is implemented both 
with and without a log transform in order to explore 
the behavior of an anamorphosis-based approach for 
maintaining positivity. Section 4 introduces our con- 
strained ensemble quadratic programming algorithm 
and section 5 shows how this algorithm performs on the 
numerical experiments introduced in section 3. Section 6 
discusses benefits and drawbacks of the proposed algo- 
rithm and identifies some open research issues. 

2. Problem formulation 

Consider a scalar quantity w whose evolution is gov- 
erned by the continuity equation with no sources or 
sinks: 

w t + V • (\w) = 0, (1) 

w(x, f Q ) = w 0 (x), for x in D, (2) 

where v is a given velocity field, to is the initial time, and 
D is the spatial domain, assumed either cyclic or to have 
no mass flux through the boundaries. Then the total 
integral of w over D is conserved through time: 

w(x,t)dx = wJx)dx , (3) 

,d )d 

and if ( ) denotes expectation, then 

(w(x,t))dx=\ (wJx))dx. (4) 

. D JD 

Here the value of the initial state w 0 at any given loca- 
tion is random but we suppose that the total mass M 
of the initial state is fixed and deterministic, so that 
Jd^oOO) dx = J D wo(x) dx = M. Note that the spatial 
distribution of mass at times after the initial time will be 
random as a result of initial condition uncertainty. We 
assume that we have access to a numerical model that 


exactly conserves a discrete version of the total mass 
integral. 

In ensemble data assimilation we typically work with 
an ensemble of spatially and temporally discretized state 
n-ve ctors that approximate realizations of w at the n 
grid points of a computational grid that covers the do- 
main D, evaluated at the time t k . The ensemble mem- 
bers are updated to incorporate information from 
observations collected at specified measurement times. 
In sequential assimilation algorithms, such as the en- 
semble Kalman filter, it is convenient to separate the 
assimilation process into two steps, carried out at each 
measurement time: 1) a forecast step that uses a nu- 
merical model of Eq. (1) to propagate the analysis en- 
semble from the previous measurement time forward to 
the current measurement time and 2) an analysis step 
that computes a new analysis ensemble from the forecast 
ensemble and current measurements. The first forecast 
in this recursion is initialized with a set of specified 
random initial conditions and all subsequent forecasts 
are initialized with the most recent analysis ensemble. 

The specific operations used to derive the analysis en- 
semble at each measurement time depend on the par- 
ticular updating procedure selected. Here we consider 
the EnKF described in Evensen (2009) and Burgers et al. 
(1998). In the EnKF update step the analysis ensemble 
member w k l is obtained by combining the forecast (prior) 
ensemble member with an m k vector of perturbed 
measurements w k \ as described by the following update 
equation: 

<■' = "L + K *«’' - ** - H rf) . ( 5 ) 

where i = 1, ... , N ens , N ens is the number of ensemble 
members, and \ k is a known possibly nonzero mea- 
surement error mean. Following usual EnKF practice, 
each perturbed measurement vector w° k ’ 1 is a random 
sample from a specified multivariate normal probability 
distribution with a mean equal to the m k vector w£ of 
actual measurements and a covariance given by the 
specified m k X m k observation error covariance matrix 
Ryt- The gain is given by 

K k=H H k(» k H H k + ^ k r\ ( 6 ) 

where is an n X m k observation matrix and P{ is the 

f K 

n X n forecast error covariance of w k . The mean 
= lW e ns2i=r wit* th e analysis ensemble, called the 
analysis, is often selected as an estimate of the uncertain 
state at t k . When the new analysis ensemble at t k is 
computed one cycle of the filter recursion is completed 
and the process repeats at t k+1 . 



February 2014 


JANJIC ET AL. 


759 


In the EnKF the forecast error covariance appearing 
in Eq. (6) is calculated as follows: 

N 

pf k = N 1 _1 X [ w fc ' - w{]K’' - w{] T , (7) 

7V ens 1 i= 1 

where w[’ 1 are the individual forecast ensemble mem- 
bers for i = 1, ... , N e ns , and is the forecast ensemble 
mean (i.e., w{ = l/N ens 

Other versions of the ensemble Kalman filter update 
include square root filters, which generate the analysis 
ensemble by first calculating the analysis ensemble mean 
and then adding a random deviation for each replicate. 
These filters do not require the use of perturbed mea- 
surements. An example is the ensemble transform Kalman 
filter (ETKF; Bishop et al. 2001; Wang et al. 2004, 2007a; 
Hunt et al. 2007). Here we use the classical perturbed 
measurement EnKF for comparison with the ensemble 
quadratic programming algorithm introduced in section 
4. This is convenient because the EnKF can be viewed as 
an important special case of the quadratic programming 
algorithm. However, the discussion of mass conserva- 
tion and nonnegativity that follows applies to all com- 
mon versions of the ensemble Kalman filter, including 
both the EnKF and the ETKF. 

From Eq. (4), we require that the continuous version, 
w a (x , tk) = (w(x, t) | wf , . . . , w£) of an analysis must 
satisfy 

f w a (x,t k )dx = f (w 0 (x))dx, (8) 

JD JD 

and that the analysis error covariance function, 

p a (x v x 2 ,t k ) 

- (M *1 ,t k ) - w a (x v t k )][w(x 2 ,t k ) - W a (x 2 ,t k )]), 

( 9 ) 

must satisfy 

P a (x 1 ,x 2 ,t k )dx 1 = 0, for all x 2 . (10) 

• D 

The latter condition reflects the requirement that every 
realization of w a (x, t k ) must conserve total mass. Such an 
analysis error covariance is called “mass conserving.” In 
the discrete case, the analysis error covariance matrix P£ 
is mass conserving if 

P*e = 0, (11) 

where e = e nX i = [11. . .1] T . The form of e given here is 
chosen for simplicity. The exact form of the definition in 


Eq. (11) will depend on the grid of our numerical model 
and the quadrature chosen for Eq. (10). 

We can define the total mass of each member in the 
forecast and analysis ensembles as M = e T w [ ,l and 
M £’ 1 = respectively. In appendix A we show that 

both the EnKF and ETKF algorithms produce mass 
conserving covariances that give the same total mass M 
for each ensemble member through the update. That is, 
if e T w^ = M for each forecast ensemble member, then 
c T w [ ,l = e T w{ — e T w^ M = e T w^ = M. 

Now suppose that w is a nonnegative scalar quantity 
such as humidity or the concentration of a chemical 
constituent. There is no guarantee that the analysis 
mean or a given analysis replicate produced by an en- 
semble filter will be nonnegative everywhere, even when 
the forecast is nonnegative everywhere and total mass is 
conserved. In fact, ensemble filters often conserve mass 
by canceling large positive values with negative values. 
This is easily shown with examples such as those de- 
scribed in the next section. Various methods, such as 
truncation of negative values to zero or formulating the 
update step in terms of log(w) can force nonnegativity 
but the resulting analysis mean and replicates typically 
no longer conserve mass. Thus, it is fairly easy to obtain 
either mass-conservative or nonnegative analyses but 
much more difficult to obtain analyses that are both 
mass conservative and nonnegative. The problem is il- 
lustrated by example in section 3 and a possible solution 
is presented in section 4. 

3. Ensemble Kalman filter performance for two 
examples 

We now consider two examples chosen to illustrate 
how mass conservation and/or nonnegativity problems 
can arise during the filter update step. These examples 
focus on feature-oriented problems where ensemble 
Kalman filters may have difficulty generating results that 
conserve mass and/or maintain the proper sign. 

a. One- dimensional static analysis with non- Gaussian 
background and observation errors 

In our first example the true feature is a static one- 
dimensional hat (isosceles triangle) function (cf. with 
the smooth function used in Lawson and Hansen 2005) 
of unit height and five grid points wide on a ID periodic 
domain. The state vector that describes this feature 
consists of values defined at 40 uniformly spaced grid 
points (see Fig. 1). The peak of the true feature is lo- 
cated at grid point 20. In a virtual experiment we 
generate noisy synthetic observations of the true state 
at a single time by summing the true (nonnegative) 
values at specified measurement locations and additive 
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Fig. 1. ID static forecast/prior. The true state (black), forecast/prior mean (red), and forecast/ 
prior ensemble (gray) are shown. Mass of true, forecast mean, and each forecast ensemble 
member is equal to 2. 


lognormal random measurement errors (also nonnega- 
tive), as follows: 

W Z = Hw fc + r ° k , (12) 

where H = is a time-invariant measurement matrix 
that consists of appropriately located zeros and ones. 
The lognormal measurement error rg has a specified 
mean (equal to 0.02 for each element) in this exper- 
iment and a diagonal covariance (variance values are 
equal to 0.01). The measurement error mean is included 
in the filter update expression, as indicated in Eq. (5). 
The spatial configuration of the measurements that de- 
termines H is discussed below. 

We assume that the exact position of the true triangle 
peak is unknown. This uncertainty is reflected through 
differences in the state vectors used to define the 50 
forecast (or prior) ensemble members, each being iden- 
tical to the true state except that the peak is located 
randomly, accordingly to a discrete uniform distribution, 
over the grid points between locations 10 and 30. Figure 1 
shows that the ensemble mean (red line) obtained by 
averaging over the forecast replicates (light gray lines) is 
nonnegative but does not preserve the triangular shape of 
the true feature (black line). 

Since neither the forecast nor measurement error dis- 
tributions are normally distributed in this example the 
Kalman filter analysis is not the exact mean of the a pos- 
teriori probability distribution and the analysis replicates 
are not samples from this distribution. Consequently, 
from a Bayesian perspective, the ensemble Kalman filter 
is suboptimal for this problem (Simon and Bertino 2009; 
Bocquet et al. 2010). Ensemble estimators are frequently 
suboptimal in applications, especially those involving 


nonnegative features with sharp boundaries. Our pri- 
mary interest here is in the filter’s ability to conserve 
mass, maintain nonnegativity, and capture the shape of 
the true triangular feature. 

With the problem setup described above we first de- 
rive the analysis replicates using the traditional EnKF 
described in Eq. (5). Figure 2 shows the true feature 
(black line) and analysis mean (red line) and ensemble 
(light gray lines) generated by the EnKF when mea- 
surements (green circles) are taken at every other grid 
point in the interval between locations 10 and 30. The 
analysis and the ensemble replicates exhibit spurious 
positive and negative lobes away from the true peak, 
although the total mass is conserved by both replicates 
and mean. The analysis ensemble replicates are gener- 
ally not positive isosceles triangles. The relatively large 
anomalies shown in Fig. 2 appear to result from poor 
interpolation of values to unmeasured locations. 

This effect is revealed in a different form in Fig. 3, 
where measurements are taken at every location over 
the more limited range from 15 to 25. In this clustered 
measurement case poor extrapolation yields large anom- 
alies outside the measured region. Figure 4 compares 
the standard deviation of the EnKF analysis ensemble 
to the root-mean-square error (RMSE) between the 
analysis replicates and the true state, for the case with 
measurement gaps. This comparison indicates that the 
analysis ensemble generally captures the true degree of 
variability, with some underestimation at unmeasured 
locations. 

The problem of maintaining nonnegativity in a Kal- 
man filter update lies in the underlying Gaussian as- 
sumptions in Eqs. (5) and (A3) (Simon and Bertino 
2009; Bocquet et al. 2010). These assumptions are not 
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Fig. 2. ID static analysis results for the EnKF with measurement gaps. The true state (black), 
observations (green), analysis ensemble (gray), and ensemble mean (red) are shown. Mass is 
conserved by all analysis ensemble members and analysis mean, but there are significant 
negative anomalies. 


applicable for estimation of a state that has discontinu- 
ities, or for a state that is nonnegative, or for problems 
characterized by an error in the location of a disturbance 
(Chen and Snyder 2007). There are several techniques 
for modifying filters to enforce nonnegativity (Cohn 1997; 
Lauvernet et al. 2009; Bocquet et al. 2010; Schneider 
1984). One alternative is to use Gaussian anamorphosis 
(Simon and Bertino 2009), which introduces a nonlinear 
change of state variables, such as a log transformation, 
in order to make the analysis step more consistent with 
Gaussian assumptions. 

We now consider a version of the EnKF, formulated 
in terms of the transformed state w k j = log(w& ; + e), 


where e is a small positive number that ensures a finite 
w kj value when w kj - = 0 and the subscript j refers to the 
;th scalar component of or w^. This transformation is 
an example of the anamorphosis approach taken by 
Simon and Bertino (2009). In our experiments we take 
e = 10 -3 . The log transformed EnKF works with a vector 
w£ of transformed synthetic measurements with com- 
ponents w£ ; = log(w£- + e). The measurement equation 
of the log transformed EnKF assumes that w£ is related 
to the log transformed state according to the fol- 
lowing additive measurement equation: 

*'k =H *'k +i °k’ ( 13 ) 



Fig. 3. ID static analysis results for the EnKF with clustered measurements. The true state 
(black), observations (green), analysis ensemble (gray), and ensemble mean (red) are shown. 
Mass is conserved by all analysis ensemble members and the analysis mean, but there are 
significant negative anomalies. 
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Fig. 4. ID static analysis results for the EnKF with measurement gaps. Comparison of en- 
semble standard deviation (solid) and RMSE between analysis ensemble members and true 
(circles). The EnKF variances are generally comparable to the RMSE, with some un- 
derestimation at scattered locations. 


where x° k is a vector of additive measurement errors in 
the observed components of the log transformed vari- 
able Wfc. This additive error equation for log trans- 
formed variables is not equivalent to the additive error 
equation for untransformed variables given in Eq. (12) 
but is required by the additive error assumption of 
the EnKF if the states and measurements are expressed 
as log transformed variables. Consequently, Eq. (13) 
should be viewed as an alternative to Eq. (12). This al- 
ternative is similar to the measurement error model 
described by Simon and Bertino (2009). 

The individual elements of x° k are assumed to be mu- 
tually independent with a uniform mean of zero and 


a uniform variance that can be adjusted to capture the 
aggregate effects of measurement error on the log 
transformed measurements (Simon and Bertino 2009). 
The need to adjust these measurement error statistics 
can be viewed as a limitation of the log transform ap- 
proach since it is difficult to know in advance how they 
should be selected. 

Results for the log transformed EnKF with measure- 
ment gaps are shown in Figs. 5 and 6. The log trans- 
formed filter is very sensitive to the variance specified 
for its assumed additive measurement error x° k . In Fig. 5 
we have set the variance of r° k equal to the variance of 
the log of x k . This value is 1.81 for the experiments 



Fig. 5. ID static analysis results for the log transformed EnKF with measurement gaps and 
additive log transformed measurement error variance set at the log of the additive un- 
transformed synthetic measurement error. The true state (black), observations (green), anal- 
ysis ensemble (gray), and ensemble mean (red) are shown. The analysis mean mass is 0.348 < 2. 
The mass is not conserved, but the analysis mean and all analysis replicates are nonnegative. 
The analysis does not capture the triangular feature. 
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Fig. 6. ID static analysis results for the log transformed EnKF with measurement gaps and 
measurement error covariance diagonal set at 0.1 variance of the log measurement error. The 
true state (black), observations (green), analysis ensemble (gray), and ensemble mean (red) are 
shown. The analysis mean mass is equal to 1.42 < 2 is not conserved, but the analysis mean and 
all analysis ensemble members are nonnegative. The ensemble is greatly constrained by the 
measurements and ensemble variability is small. 


described here. The variance used to generate Fig. 6 is 
0.181, one-tenth of the value used in Fig. 5. The analysis 
values in the higher measurement error variance case of 
Fig. 5 are small and the measurements have relatively 
little impact. The mass of the analysis is 0.348, which 
is much lower than the true mass of 2. The lower mea- 
surement error variance case shown in Fig. 6 is highly 
constrained by the measurements, with minimal en- 
semble spread and an analysis mass of 1.42. This value is 
still significantly lower than the true mass. 

In the limit of very low measurement error variance 
the log transformed EnKF analysis ignores the fore- 
cast, giving a triangle with a smaller base and smaller 
total mass than the true. The mass conservation results 
are poorest for intermediate values of the measure- 
ment error variance, when the analysis deviates from 
the forecast but fails to capture the shape of the true 
feature. Similar results are obtained for the clustered 
measurement configuration, which is not shown here. 
When the measurement error random seeds are varied 
some of the analysis ensemble members can take on 
positive values significantly higher than 1.0, even for 
small measurement error variances. Generally speak- 
ing, the log transform EnKF appears to be marginally 
stable for this problem, requiring fine adjustments of 
the measurement error variance to give reasonable 
analysis results. In any case, although analysis ensem- 
ble members and the analysis are nonnegative every- 
where, mass is generally not conserved. 


The results for this simple one-dimensional example 
show the dilemma encountered with classical ensemble 
Kalman filtering for problems where mass conservation 
and sign are both important. The classical EnKF con- 
serves mass if the forecasts are mass conservative but it 
can produce unrealistic analyses (mean and ensemble) 
that are negative. The log transformed filter gives non- 
negative analyses but does not generally conserve mass. 

b. Two-dimensional dynamic analysis 

The previous example shows the behavior of the 
EnKF solution for only one measurement update. To 
estimate the effect of analysis errors through time and 
on forecasts, we perform a second virtual experiment that 
considers two-dimensional solid-body rotation (Tremback 
et al. 1987; Janjic et al. 2011). In this experiment a moving 
cone completes a full clockwise rotation of 2i t about the 
origin (domain center) every 48 h. Synthetic observa- 
tions are assimilated every quarter revolution (4320 time 
steps) until seven forecast/analysis cycles are completed 
(time step 30241). The experiment is carried out on 
a uniform numerical grid of 101 by 101 square cells, each 
of size 8 km by 8 km. 

The initial ensemble is specified to be a set of cones, 
each with a radius of 100 km at the base and a height of 
100 units. The true initial cone is centered at the grid 
point with indices i = 33 and j = 33. The central location 
of each cone in the 50-member ensemble is described by 
an angle and radius defined relative to the domain 
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Fig. 7. Three members of the initial 50-member ensemble generated by perturbing the location 

of the center of the cone. 


origin. The angle is perturbed around the true cone 
angle by a random number normally distributed with 
mean 0 and standard deviation of 0.5 radians. The dis- 
tance of the cone from the center of rotation is perturbed 
as well, with normally distributed random error with 
mean 0 and standard deviation of 0.5 km. The random 
angle and the distance define the center of each cone in 
the ensemble. The ensemble of cones is transported with 
pure advection by simply changing the angle of their 
centers relative to the domain origin at the appropriate 
rate. This ensures that the solution to the governing 
equation given in Eq. (1) is exact, with no numerical 
dispersion or other numerical errors. Three of the en- 
semble members are shown at the initial time in Fig. 7. 

The synthetic observations are obtained by perturbing 
the true solution with a small amount of log normally 
distributed noise with mean 0.5 and variance 1. The 
minimum observation value is 0.002. In this experiment 
the observation operator varies in time. It uses every 
20th synthetic observation at the locations where true 
cone would have values greater than zero. This corre- 
sponds to measurements at 25 grid points per analysis 
time. The restriction of noisy measurements to the re- 
gion of the true cone reflects the fact that measurements 
outside this region will be small compared to the cone 
amplitude and would normally be removed by truncat- 
ing measurements below an appropriate threshold. 

The forecast mean at the first update time is computed 
as an average over the ensemble members, producing 
a much attenuated field with a maximum value of 46.7 
and minimum value of zero (Fig. 8). Note that, although 
every ensemble member generated at the initial time is 
a perfect cone with the desired properties, the structure 
of the cone is not preserved in the ensemble average. 
Since the model is exact, each copy of the cone rotates 
without losing its structure or its minimum and maxi- 
mum values. However, the shape of the cone is not 


preserved through the EnKF analysis step. Although the 
total mass is conserved, unrealistic negative mass values 
are obtained after the analysis. Figure 9 shows the 
analysis and the three analysis ensemble members at 
the end of the experiment (time step 30 241) from bird’s- 
eye perspective. The minimum, maximum, and RMS error 
values of the analysis are —11.2, 98.4, and 1.1, respectively. 
The analysis and each of the ensemble members show 
spurious positive and unphysical negative values away 
from the cone structure. Negative values (depicted with 
dark blue contour lines in Fig. 9) reach - 16.3 in the en- 
semble. The ensemble members that have the lowest 
(—16.3) and the highest (—7.3) minimum values are shown 
in Fig. 9, together with an ensemble member with mini- 
mum value of -10.5. The difference between the chosen 
ensemble members is the largest in the area away from 



Fig. 8. The forecast ensemble mean at the first measurement time 
from bird’s-eye perspective. 
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Fig. 9. (top left) The analysis at the end of the solid-body rotation experiment (time step 30241), obtained with 
the EnKF algorithm and 50 ensemble members. Ensemble members are shown as examples of replicates with the 
(top right) highest and (bottom left) lowest minimum values, (bottom right) An ensemble member with a minimum 
value between the two is depicted. Contour lines in the range from -10 to 10 are shown in steps of 1, and above 10 in 
steps of 10. 


the cone, where spurious error structures appear (see 
Fig. 9). The maximum value of the true cone (100) is 
slightly underestimated in the ensemble since the max- 
imum of the ensemble members varies between 96.9 and 
100.5. Although all ensemble members show the cone 
structure in the right location, unphysical negative values 
away from the cone are large. Both positive and negative 
spurious values affect a large area of the domain. 

4. Problem solution 

The update step of the ensemble Kalman filter de- 
scribed in Burgers et al. (1998) can be posed as the so- 
lution to a set of regularized least squares optimization 
problems, one problem for each member in the ensem- 
ble. The problem associated with members i = 1, , 

N e ns at time t k can be expressed using the same notation 
as in Eqs. (5) and (6): 


w ^ = w ^ + ar gmin^[5w iT (P / ) 1 5w i + f T R 1 f l ], 

8w l ^ 

(14) 

where — w{’* is the analysis increment, which 

serves as the decision variable, and f = - 

r k = w- — H&w{ ,z — - r£. Note that is obtained 

from the forecast ensemble, as indicated in Eq. (7). Also, 
the measurements appearing in f* are perturbed as in Eq. 
(5). The solutions to the A ens optimization problems de- 
fined in Eq. (14) form the analysis ensemble. When there 
are no constraints in Eqs. (5) and (6) give a closed form 
solution to the ensemble optimization problem in Eq. (14) 
(Zupanski 2005; Wang et al. 2007a). 

Our extension of this optimization formulation of 
the ensemble Kalman filter adds linear inequality con- 
straints in order to enforce nonnegativity in the update. 
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The resulting constrained analysis step solves the fol- 
lowing quadratic programming problem for each mem- 
ber i — 1, . . . , N ens . 

w"-' = w(-' + argmin 1 [5w' T (P / ) _ 1 5w'' + f'R 'f] 

(15) 

subject to the following nonnegativity constraint: 

Sw'>-wf’\ (16) 

Note that the inequality in Eq. (16) is equivalent to 
w£ l > 0. In general the constrained solution to Eq. (15) 
cannot be expressed in closed form but must be derived 
numerically. 

The quadratic programming objective given above 
depends on the inverse of the sample covariance P^, 
which is singular when N ens — 1 < n and/or the forecast 
covariance is mass conserving. The low rank of P ; al- 
lows us to transform the optimization problem by re- 
ducing the number of decision variables to p = Rank(P^), 
which is no larger than N ens — 1. Appendix B shows that 
the quadratic programming solution conserves mass if 
the forecast covariance is mass conservative and the 
decision variable 5w* is chosen to lie in the p di- 
mensional subspace spanned by the forecast ensemble. 
For the unconstrained case this solution duplicates the 
closed form EnKF solution given in Eq. (5). 

If 8w l lies in the forecast ensemble subspace it can be 
expressed in terms of a p-dimensional transformed de- 
cision variable r\ as follows: 

8w l = Lr] 1 , (17) 

where the columns of the n X p dimensional matrix 
L form a basis for the forecast ensemble subspace and 
the elements of rf can be viewed as the weights in 
a linear combination of the basis vectors. Then P^can 
be written in terms of the p X p dimensional covari- 
ance Q of 77 : 

P / = LQL t . (18) 

Since there is flexibility in defining L we select it to 
satisfy the requirement that P^ = LL T , so that Q is the 
p-dimensional identity matrix and L is the matrix square 
root of P^. This square root may be computed using 
a singular value decomposition of the matrix whose 
columns are the differences between the vectors of the 


forecast ensemble members and the forecast ensemble 
mean. 

We use the change of variables given in Eq. (17) to 
rewrite the constrained quadratic programming prob- 
lem in terms of r\\ 

i? ! = argmin firj'V + f ,T R _ 1 f']> (19) 

V 2 

where f = w£’* - - H k Liq l - x° k subject to the fol- 

lowing nonnegativity constraint: 

( 20 ) 

The analysis ensemble can be derived from the solution 
to Eq. (19): 

w k ‘ = w i f k + Sw ' = w k‘ + w • (21) 

The analysis ensemble and its mean all lie in the forecast 
ensemble subspace since w[ J and 8W lie in this space by 
construction (see appendix B). 

Once the analysis ensemble is calculated using 
Eq. (21), it is propagated with the forecast model to 
obtain the new forecast ensemble, as in the classical 
unconstrained ensemble Kalman filter. For easy refer- 
ence, we call this ensemble data assimilation method 
QPEns. The structure of the QPEns algorithm insures 
that the analysis and each member of the analysis en- 
semble will have the desired mass conservation and 
positivity properties if the forecast ensemble is mass 
conservative. 

There are a number of QP algorithms available to 
compute the optimum for Eqs. (19)-(20). For example, 
the active-set method can be viewed as an extension of 
the traditional EnKF analysis. The first iteration starts 
with an unconstrained optimization. If this solution 
satisfies all constraints, then the optimum is found. 
If not, more iterations will follow, where in each it- 
eration some inequality constraints are converted to 
equality constraints (i.e., made active) or removed as 
equality constraints (i.e., made inactive). The number 
of iterations depends on the problem at hand. For 
the experiments in section 5a, we found that 10 iter- 
ations are sufficient on average. For many applications, 
obtaining the ensemble forecast with the forecast model 
will dominate the computations, so the additional 
effort required by the QPEns optimization will likely 
be affordable. There is much potential in future re- 
search for developing more efficient optimization pro- 
cedures that exploit the special structure of the QPEns 
problem. 
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Fig. 10. ID static analysis results for QP with positivity constraint and nominal assumed 
measurement error variance 0.01. True state (black), observations (green), analysis ensemble 
(gray), and analysis ensemble mean (red). The analysis ensemble members and their mean 
conserve mass and are always nonnegative. The analysis mean underestimates true at the peak 
in order to conserve total mass since the ensemble members and mean both include some mass 
in the region where the true state is zero. 


5. Quadratic programming performance for two 
examples 

a. One- dimensional static analysis with non- Gaussian 
background and observation errors 

This section considers the static one-dimensional 
problem discussed in section 3 when the analysis en- 
semble is derived with the ensemble quadratic pro- 
gramming algorithm (QPEns). All results were obtained 
with an active set Matlab quadratic programming routine. 
Figure 10 shows the true state (black) and the analysis 
mean (red line) and ensemble (light gray lines) generated 
by QPEns for the measurements already considered in 
Fig. 2 (green circles). In this case the measurement error 
variance used in the QPEns is set equal to the variance of 
the additive lognormal synthetic measurements, without 
any adjustment. For our computational experiment the 
mean and variance of the lognormal measurement errors 
are 0.02 and 0.01, respectively. These are the same values 
used to generate the untransformed EnKF results pre- 
sented in section 3. The resulting QPEns analysis mean 
and ensemble all conserve mass and are nonnegative ev- 
erywhere. Some of ensemble members take on nonzero 
values in regions where the true field is zero, leading to 
nonzero analysis values in these regions. As a result, the 
analysis mass near the true peak needs to decrease in 
order to maintain the correct total analysis mass (i.e., the 
areas under the red and black curves need to be the same). 

The QPEns results are sensitive to the specified mea- 
surement error covariance R. Figure 11 shows the results 


obtained for the same problem considered in Fig. 10 
with the measurement error variance used in the QPEns 
algorithm reduced to one-quarter of the actual value 
(i.e., from 0.01 to 0.0025). In this case, the measurements 
have more influence and the analysis mean and all the 
ensemble members are very close to the true values. 
This behavior reflects the fact that, for this particular 
problem, the low measurement error covariance qua- 
dratic programming objectives are minimized when the 
values at the peak and the two nearest measurement 
locations are as close as possible to the corresponding 
measurements. Mass conservation and nonnegativity 
can only be satisfied if the values at points farther from 
the peak are all close to zero. It is interesting that the 
EnKF still gives significant negative values (not shown 
here) for the low measurement error variance case. In 
this example the nonnegativity constraint that distin- 
guishes the QPEns algorithm has an important impact 
on overall accuracy as well as sign. 

Figure 12 compares the QPEns standard deviation of 
the analysis ensemble to the RMSE between the analysis 
ensemble and the true state, for the nominal specified 
measurement error variance value. The analysis ensem- 
ble variances are generally comparable to the RMSE, 
with some underestimation at a few locations near the 
true peak. Note that the RMSE values for QPEns are 
significantly lower than the EnKF values plotted in Fig. 4, 
indicating a better match to the true feature. 

The results shown above indicate that the constrained 
QPEns algorithm provides a major improvement in 
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Fig. 11. ID static analysis results for QPEns with positivity constraint and lowered assumed 
measurement error variance (0.0025). True state (black), observations (green), analysis en- 
semble (gray), and analysis ensemble mean (red). The analysis ensemble and their mean 
conserve mass and are always nonnegative. The analysis mean and ensemble members are very 
close to the true values. 


accuracy and physical realism over either an un- 
transformed or log transformed EnKF, at least for the 
problem we have considered. The QPEns conserves 
mass, gives nonnegative values, and provides an accu- 
rate description of the true feature of interest. 

b. Two-dimensional dynamic analysis 

In this section we use the QPEns algorithm to solve 
the two-dimensional dynamic data assimilation problem 
introduced in section 3b. For implementation of the 
QPEns algorithm in this example, we constructed the 
perturbed observations as w£ + r£’\ where r] V is a vector 
normally distributed A/*(0, R&) and Rk is a diagonal 
matrix with l 2 on the diagonal. The results of this 


experiment are summarized in Fig. 13. The QPEns al- 
gorithm is able to recover the cone structure, with 
a maximum value of 94.9 and RMSE of 0.4 at the end of 
the experiment (cf. Fig. 9). As in Fig. 9 we show an ex- 
ample of three ensemble members. Since the result of 
the QPEns cannot be negative, we show the ensemble 
members with lowest and highest maximum values, and 
one with the maximum value in between. Maximum 
values of ensemble members vary in the range between 
56.6 and 98.5. One of the depicted ensemble member 
(Fig. 13, bottom-left panel) almost perfectly represents 
the true cone structure with the maximum value of 98.5. 
The ensemble member with the lowest maximum value 
differs from the true cone, with the errors primarily 



Fig. 12. ID static analysis results for the QPEns with positivity constraint and nominal 
measurement error variance. Comparison of ensemble standard deviation (solid) and RMSE 
between analysis ensemble members and true (circles). The QPEns variances are generally 
comparable to the RMSE, with some underestimation at a few locations near the true peak. 
RMSE values for the QPEns are significantly lower than for the EnKF (cf. Fig. 4). 
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Fig. 13. (top left) The analysis at the end of the solid-body rotation experiment (time step 30241), obtained with 
the QPEns algorithm. Ensemble members are shown as examples of replicates with the (top right) lowest and 
(bottom left) highest maximum values, (bottom right) An ensemble member with maximum value between the two is 
depicted. Contour lines in the range from -10 to 10 are shown in steps of 1, and above 10 in steps of 10. 


localized around the cone in the range from 0 to 10. The 
third ensemble member shown in Fig. 13 has a maximum 
value in between the lowest and highest. This member 
has a clearer cone structure than the one with the lowest 
maximal value but positive errors are still present around 
the cone. In all cases, the analysis ensemble members 
have the same mass as the true cone and are nonnegative, 
resulting in a more accurate analysis mean than the tra- 
ditional EnKF. 

6. Discussion and conclusions 

In this paper we propose using constraints to enforce 
mass conservation, nonnegativity, and other physical re- 
quirements in an ensemble Kalman filter update. Fore- 
cast means and covariances convey useful but sometimes 
incomplete information about the physical requirements 
imposed by conservation laws. As a result, ensemble 
Kalman filter updates that rely only on these moments 
and scattered noisy measurements can produce unphysical 
analyses and analysis ensemble members. It is possible 


to deal with some physical requirements, such as mass 
conservation, through proper construction of the 
forecast error covariance and the subsequent update. 
But this does not generally insure that analysis results 
are nonnegative. Conversely, it is possible to impose 
nonnegativity by using transform methods such as 
anamorphosis, but these methods do not necessarily 
conserve mass. If the classical unconstrained ensem- 
ble Kalman filter update is appropriately constrained 
it is possible to conserve mass and also to maintain the 
correct sign. 

When measurements are linearly related to the state, 
the ensemble Kalman filter update can be posed as a set 
of unconstrained quadratic programming problems, one 
for each replicate. The solutions to these unconstrained 
problems can be expressed in closed form. The qua- 
dratic programming structure of the problem is main- 
tained if linear equality and inequality constraints are 
added, but the solutions must generally be obtained from 
a numerical optimization procedure rather than from 
a closed form expression. Fortunately, many important 
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physical constraints, including mass conservation and 
nonnegativity, are linear and fit into a quadratic pro- 
gramming framework. This makes it convenient to add 
constraints to existing sequential data assimilation algo- 
rithms based on ensemble Kalman filters. The quadratic 
programming formulation has a number of important 
advantages, including the availability of very efficient 
solution algorithms and the guarantee that the problem 
has a unique minimum for linear observation operators 
when the Hessian of the quadratic objective function is 
positive definite. 

The benefits of including constraints to enforce non- 
negativity are apparent in the results obtained in our 
two synthetic experiments. In both cases, ensemble 
quadratic programming captures the shape and mass of 
a distinctive spatial feature through a filter update with 
noisy measurements. The quadratic programming ap- 
proach works much better than a classical ensemble 
Kalman filter, which gives negative estimates even though 
all measurements and forecasts are nonnegative. The 
classical EnKF conserves total mass by generating in- 
flated positive masses in some locations in order to 
cancel negative masses generated in other locations. The 
log transformed EnKF is able to maintain nonnegativity 
but does not generally conserve mass and requires ad- 
justment of its measurement error covariance in order to 
obtain reasonable ensemble members and to capture the 
approximate shape of the true feature. 

The quadratic programming approach described here 
has some limitations that are important to note. It is 
appropriate when the observation operator is linear and 
when all constraints included in the update are linear. 
The forecast model can be nonlinear, as in other en- 
semble filtering methods that compute forecast statistics 
by propagating ensemble with nonlinear dynamic models. 
The basic concept of constraining the Kalman filter 
update can be extended to accommodate nonlinear 
measurement operators and constraints. However, the 
optimization must then be performed with a more ex- 
pensive nonlinear programming algorithm and there is 
no longer a guarantee of a unique minimum. It is pos- 
sible that the quadratic programming formulation adop- 
ted here could be retained for nonlinear measurement 
operators if the objective function to be minimized is 
expressed in terms of the analysis error covariance rather 
than the forecast and observation error covariances 
(Zupanski 2005). Linear equality and inequality con- 
straints could then be included as described in section 4. 

As mentioned above, the ensemble quadratic pro- 
gramming approach requires numerical solution of a 
different quadratic programming problem at every 
analysis time, for every ensemble member. This is not 
a significant limitation for our simple examples but it 


could require substantially greater computational 
effort than the standard closed form ensemble Kalman 
filter update, especially for spatially distributed prob- 
lems with many degrees of freedom. However, the 
relative increase in overall computational effort may 
not be significant because in many high-dimensional 
ensemble filtering problems computational effort 
is dominated by the forecast rather than the analysis 
step. Larger problems will have to be investigated 
before we can assess the overall impact of solving a 
quadratic programming problem for every analysis 
ensemble member. It may be possible to reduce com- 
putational effort by taking advantage of the fact that 
quadratic programming solutions for the different en- 
semble members tend to be clustered around a common 
mean. Also, the different quadratic programming solu- 
tions required in the analysis step can be computed in 
parallel. 

Our quadratic programming approach for including 
constraints in an ensemble Kalman filter is related to 
other ensemble data assimilation methods including 
hybrid variational methods and randomized maximum 
likelihood. These various approaches are complemen- 
tary and it is likely that they could be combined in var- 
ious ways. The distinctive aspects of our approach are an 
emphasis on the need to include physical constraints 
during the update and a formulation that takes advantage 
of the computational benefits of quadratic programming. 
Together, these provide a practical and effective way to 
ensure that data assimilation results satisfy fundamental 
physical requirements. 
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APPENDIX A 

Mass Conservation in Ensemble Kalman Filters 

In ensemble Kalman filters the sample forecast error 
covariance matrix in Eq. (7) is derived from an ensemble 
of states produced by integration of a numerical model. 
Many numerical integration schemes can conserve the 
total (global) mass of tracers (e.g., Schneider 1984; Lin 
and Rood 1997). Each replicate in a forecast ensemble 



February 2014 


JANJIC ET AL. 


771 


computed with a conservative scheme conserves total 
mass and the sample covariance derived from the fore- 
cast ensemble is mass conserving. That is, e T w£’* = M 
for each ensemble member and the mean w£ over any 
number of ensemble members has the same mass M, so 
e T (wf - w{) = 0 for all i and Eq. (7) gives e T P{ = 0. 

Since e T P{ = 0, Eq. (6) gives e T K^ = 0 and it follows 
from Eq. (5) that e T w^ = e T w^ = M for the classical 
EnKF. Therefore, the analysis ensemble and analysis 
mean of the EnKF all conserve mass if the forecast en- 
semble members conserve mass or, equivalently, if the 
forecast covariance is mass conserving. In this case the 
EnKF analysis covariance is also mass conserving. This 
covariance can be expressed, for any K k , as 

K = (I - K,H a ,)P{(I - K A H fc ) T + k a r a k l (Al) 

Since P[ is mass conserving and the gain satisfies 
e T K k = 0, then from Eq. (Al) we have that e T P a k = 0, so 
P a k is mass conserving. 

Similar reasoning applies to the ETKF. In this case the 
analysis mean is computed directly from 

< = w ( + K 4< - ~ H k w l) • ( A2 ) 

Consequently, it follows that e T w£ = e T w{. This in- 
dicates that the ETKF analysis is mass conserving if the 
forecast model conserves mass. 

The ETKF generates its analysis ensemble as follows: 

< ,i = <+ V^ens-l [Wfrji’ A 3 ) 

where Wf = l/\AVens - I fw /' 1 - y/{, . . . , wf ,iVens - w{] is 
the n X A ens matrix of deviations of the forecast en- 
semble members from their mean. We take the A ens X 
N e ns transformation matrix as in Wang et al. (2004, 

2007a) and Hunt et al. (2007). Namely, if C k and are 
the matrices of eigenvectors and corresponding eigen- 
values of the matrix = C*D*(C fc ) T , 

respectively, then T/ { = C^(l^ ens + Dk)~ 1,2 Cj is the ma- 
trix that transforms deviations of forecast ensemble 
members from the forecast mean into deviations of 
analysis ensemble members from the analysis mean. 
Here, \ Nens denotes the A ens X A ens identity matrix. 

The ETKF analysis error covariance can be shown to 
be P a k = W/T,j[wf (Bishop et al. 2001; Wang et al. 
2007a). It follows that the ETKF analysis covariance is mass 
conserving since e T P;:=e T W[T,Tjw( and e T W{ = 0. 
Also, from Eq. (A3) we have e T w A ' = e' w)' = M . 

The above discussion indicates that the analysis en- 
semble and analysis mean of both the EnKF and the 
ETKF conserve mass and their analysis covariances are 


mass conserving if the forecast ensemble members 
conserve mass or, equivalently, if the forecast co- 
variance is mass conserving. For this derivation we made 
no assumptions on linearity of the model dynamics. The 
dynamics can be nonlinear as long as the numerical 
discretization scheme conserves mass. 

The proof that the analysis error covariance is mass 
conserving if the forecast error covariance is mass con- 
serving is applicable to any mass conserving forecast 
error covariance, not only those derived from an en- 
semble. Another example of a covariance formulation 
that would conserve mass is one obtained by eigenvalue 
decomposition on a sample that has a constant spatial 
integral, since then e would be an eigenvector corre- 
sponding to a zero eigenvalue of the covariance matrix 
computed from the sample, and, therefore, would be 
orthogonal to the other eigenvectors that are part of the 
low-rank modeled covariance. Since this is a usual tech- 
nique for initializing ensemble square root Kalman filter 
algorithms, we assume that we start initially with a mass 
conserving covariance matrix. Note that if the forecast 
error covariance is modeled in such a way that all its el- 
ements are positive [Gaussian, second order autore- 
gressive (SOAR), third order autoregressive (TOAR)] or 
nonnegative, then it cannot be mass conserving, since e 
cannot be a null vector of a matrix with all positive 
elements. Similarly, if localization is applied to the 
ensemble -derived forecast error covariance through a 
Schur multiplication, then the mass will not be conserved. 

APPENDIX B 

Mass Conservative Properties of the QPEns Analysis 

The mass conservation properties of the QPEns anal- 
ysis are related to the properties of the p-dimensional 
subspace spanned by the ensemble of A ens forecast en- 
semble members. This subspace is also spanned by the 
ensemble of forecast deviations w k l - w{, since the 
mean w[ also lies in the same subspace as the forecast 
ensemble members. The analysis ensemble produced by 
the QPEns algorithm or, equivalently, the analysis de- 
viations w k l - w k , are constrained by construction to 
also lie in the forecast ensemble subspace [see Eq. (17)]. 
Since the analysis deviations lie in the subspace spanned 
by the forecast deviations, they may be written as a lin- 
ear combination of the following form: 

N 

ens 

W fc' - < = W fc), ( B1 ) 

/= 1 

where the cqy, i — 1, , N ens and j = 1, , N ens are the 
scalar coefficients of the linear combination for replicate 
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i. This implies that the analysis ensemble is mass con- 
servative if the forecast ensemble is mass conservative 
since 


N 


e T ( w ^ — W D = e T X a ij( yv k ~ w {) 

7=1 


N 


= X ^ e T (w£’ 7 - w£) = 0. (B2) 

7=1 


The final equality is a mathematical statement of the 
assumption that the forecast ensemble is mass conser- 
vative. Note that this result applies to any quadratic 
programming problem with a vector decision variable 
that is constructed to lie in the forecast ensemble sub- 
space, as is done in Eq. (17). 
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